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1. SUMMARY 


The magnetic lineations discovered by MGS have been considered to be evidence of early plate 
tectonics on Mars. However, the lineations approximately follow lines of latitude, i.e., small 
circles. This presents significant geometrical problems for plate-like spreading, particularly at 
high latitudes. However, the sublatitudinal orientation of the lineations is consistent with 
meridianal extension and perhaps limited crustal spreading due to a stress event centered near the 
geographic pole. We hypothesize that this event was the early formation of the crustal 
dichotomy through mantle-convective processes. This could have taken the form of a southern 
megaplume that formed the thick highlands crust or as subduction or downwelling in the north. 
Both would result in tensional stresses in the south that would form extensional fractures 
perpendicular to the CM-CF offset. The observed magnitude and distribution of magnetization 
indicates that crustal intrusion associated with this major mantle-convective event resulted in 
-1000 km of extension in the southern highlands. Subsequent spin-axis reorientation due to loss 
of crust in the north or gain of crust in the south brough the CM-CF offset into its present N-S 
alignment. A portion of the ancient valley networks observed in the southern highlands are 
spatially associated with crustal magnetism and are quantitatively shown to be consistent with 
hydrothermal discharge over crustal intrusions. 
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3. INTRODUCTION 


The remarkable magnetic anomalies discovered by the Mars Global Surveyor {Acuna et al., 
1999; Figure 1) strongly suggest that early Mars had an intrinsic magnetic field and that newly 
created crust was magnetized by thermoremanence. The large magnitude and size of the 
anomalies further indicates that the tectonic processes that influenced their formation were global 
in scale. Connerney et al. (1999) have suggested that this process was plate tectonics, a 
hypothesis that has received widespread attention. They based this hypothesis qualitatively on 
the apparent linear nature of the magnetic anomalies, which was considered to be similar to 
terrestrial magnetic "stripes" formed by seafloor spreading. Connerney and co-workers further 
showed through simple modeling that the intensity of magnetization had to be relatively strong 
compared to terrestrial analogs in order to reproduce the large magnetic fields observed from 
orbit. 



-1500 0 1500 

Fieure 1. Magnetic field measured by MGS, from Connerney et al. (1999). Note strong sublatitudinal orientation to 
the lineations. Solid lines and dots denote segments used for geometric analysis (see below). 
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Plate tectonics has been such a successful theory on Earth because, like all good theories, it 
has strong predictive power. These predictions follow from the simplicity of rigid-plate motion 
on a sphere, which constrains the geometry of plate interactions. In general, spreading centers 
are not constrained to follow any particular orientation, although a great-circle path is required 
for normal spreading (when movement is perpendicular to the segment orientation). A second 
basic observation that can be made in addition to their linearity is that the anomalies 
approximately follow lines of latitude (Figure 1.), i.e, small circles. Especially at high latitudes, 
this presents strong geometrical problems for plate spreading, and suggests that other, alternative 
hypotheses should be considered. 

Any such alternative process must be able to produce dominantly meridianal extensional 
stresses over at least the latitudinal band in which the magnetic anomalies are observed, large 
strain (presumably tens to hundreds of percent or more; see below), and be limited to the 
Noachian era (as the lineations are observed only in the heavily cratered southern hemisphere). 
Potential candidates are global expansion, true polar wander, despinning, and formation of the 
"hemispheric dichotomy." While expansion produces extensional stress, such stress is not 
limited to a meridianal direction; furthermore the strain for plausible expansion is very small (see 
Solomon, 1979). Polar wander can be rejected for similar reasons (see Melosh, 1980; Grimm 
and Solomon, 1 986). Despinning does produce a zone of east-west normal faulting at latitudes 
higher than 36° ( Pechmann and Melosh, 1979), but the martian anomalies exceed this range, and 
plausible changes in flattening still produce total hemispheric extension limited to tens of 
kilometers. 


Mars is divided by a "hemispheric dichotomy" into heavily cratered southern highlands and 
resurfaced northern lowlands. Leading hypotheses for this dichotomy encompass both internal 
(mantle convection: Lingenfelter and Schubert, 1973; Wise et al., 1979) and external (impact: 
Wilhelms and Squyres, 1984) mechanisms. However, MGS measurements of martian 
topography (Smith et al., 1999) reveal that much of the difference in elevation between the south 
and north is due to the center-of-mass/center-of-figure (CM-CF) offset. This degree-one offset 
corresponds to a steady downward slope of 0.036° between the south and north poles, and 
strongly favors a "smooth" internal process for the formation of the hemispheric dichotomy. We 
suggest that the north-south CM-CF offset is not coincidence, but the result of spin-axis 
reorientation due to mantle-convective formation of the hemispheric dichotomy, and that the 
orientation of magnetic anomalies along lines of latitude is a result of meridianal extension of 
hundreds of kilometers or more associated with this major geological event in the early history of 
Mars. 

The objectives of this project were to (1) Quantify the planform and intensity ot crustal 
magnetization, (2) Assess how these geometrical patterns constrain models for the formation of 
the magnetization, particularly the origin of the crustal dichotomy, (3) Explore the relationships, 
if any, of crustal magnetization to surface geology', in particular, the valley networks. 
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4. MAGNETIC MODELING 

Inversion of the pattern of crustal magnetization from magnetic anomalies is the foundation 
for determining the time scale for geomagnetic reversals and measurement of seafloor spreading 
{Vine and Matthews, 1963). The crustal magnetism is preferred to the magnetic anomalies 
themselves because superposition of the anomalies of adjacent blocks with differing 
magnetization intensities and/or directions results in a complex signal that may not be visually 
interpretable; furthermore induced moments distort the observed field and the inferred positions 
of reversals. Although the effects of induced magnetism may be ruled out for Mars (its crustal^ 
magnetism must be wholly remanant in the absence of a present-day intrinsic held), reduction of 
the observed data to the actual magnetization is an important step in understanding the origin of 
the magnetic anomalies. 

Connerney et al. (1999) presented simple models for several 
individual MGS orbits. Their model is based on the following 
assumptions: (1) magnetization is two-dimensional, (2) 

magnetization is constant in crustal blocks 200 km wide by 30 
km thick, and (3) the two components of magnetization J x and J z 
are allowed to vary independently. These models reproduce the 
observed anomalies quite accurately (Figure 2). As 
acknowledged by these authors, however, this is only one 

possibility for the actual crustal magnetization 

Figure 2 . Horizontal (top) and vertical (bottom) magnetic anomalies with 
inferred horizontal and vertical crustal magnetization (middle), from 

Connerney et al ( 1 999). 

We independently modeled two of the orbits presented by 

Connerney et a/. (1999). The data were digitized from the 

published paper. We also used a generalized inverse using 
. ,«» .ooo »«o singular-value decomposition to solve for the crustal 

*<*"> magnetization. Our results for orbit 1999 Day 15 P6, also using 

an unconstrained J x J z model, are very similar to those of Connerney and co-workers (Figure 3). 
This result was obtained for a relatively low damping factor (SV cutoff ratio) of 10' . A similar 
comparison was obtained for orbit 1999 Day 20 P7. 

We performed a more thorough exploration of the parameter space for this model than 
reported by Connerney et al. (1999). We are examined the effects of varying block thickness, 
block width, and the damping factor on the goodness-of-fit, the model roughness, and the 
maximum magnetization. We found, as expected, that goodness-of-fit and model roughness are 
inversely related through the damping parameter. Because Connerney and co-workers 
emphasized the best fits to the data, their models are relatively rough. Connerney et al. (199 > ) 
state that "reversals in the direction of magnetization of adjacent blocks are common, if not the 
rule." Flowever, strong variations between adjacent model blocks is almost always a signature o 
excessively high model variance or roughness. Reversals would be more believable if groups o 
at least several model blocks had the same sign, which can be achieved by using a narrower 
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block and a higher damping. We find that goodness-of-fit is insensitive to block thickness but 
does improve with narrower blocks (because there are more free parameters) but that, conversely, 
model roughness does not vary as strongly with block width as it does with block thickness. The 
maximum crustal magnetizations, already apparently extreme for the nominal model of 
Connemey and co-workers, increase as the blocks are thinned. 

While the models of Connemey et a 1. (1991) were useful as exploratory tools to describe the 
general magnitude of crustal magnetization, the spatial patterns ol this magnetization are random 
and unrealistic (Figure 3). This is because J x and J z are both allowed to vary independently. In 
reality, the magnetic inclination will vary with absolute plate motion and true polar wander, but 
is often considered constant for some number of geomagnetic reversal cycles of seafloor 
spreading. Alternatively, if there is no plate motion and polar wander, the magnetizing field 
should closely approximate a global dipole. We have tested this latter hypothesis by constraining 
j x /j z = o.5 tan(0), where © is the geomagnetic colatitude. The magnetization amplitude (positive 
or negative) for each block remains as a free parameter. The solved magnetizations now vary 
smoothly (Figure 5) and, because the model is more tightly constrained, the required intensity of 
magnetization is less than half that of the free J*,J Z model. The cost is a modest decrease in 
goodness-of-fit (Figure 6) which, while statistically significant, does not detract from the 
plausibility of the model. 
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Figure. 3 . Best-fitting crustal magnetization 
for orbit 1999 Day 15 P6. The model 
formulation closely follows that of 
Connerney et al. (1999). The magnetization 
directions of each block are unconstrained, 
i.e., the x- and z- magnetizations vary 
independently. Note close match to the result 
of Connerney et al. (1999) (Fig- 2 above). 
However, when magnetizations are plotted as 
vectors, a random, unrealistic magnetizing 
field is inferred. 



gure 4. Data (symbols) and 2D model fit (lines) to magnetic field for model in previous figure. Goodness-of-fit is 
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Figure. 5 . Best-fitting crustal 
magnetization for same data as above, 
except magnetization directions are 
constrained to follow a global dipole. 
The latitude of the magnetic pole is a 
new free parameter, but each of the 
model blocks now has only one free 
parameter (the total magnetization) 
instead of two. The maximum 
magnetization in the constrained model 
is about half that of the unconstrained 
model (20 vs 42 A/m). The apparent 
magetic equator is near x=~1000 km 
(where the magnetizations are nearly 
horizontal) and corresponds to a 
magnetic pole near 54 °N, 16° W. 



Figure 6. Data (symbols) and 2D model fit (lines) to magnetic field for model in which magnetization direction is 
constrained to follow* a global dipole. Goodness-of-fit is 94.5 /o. 
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While the 2D models have lead to useful understanding of the crustal magnetization of Mars 
(see below), continuity between orbit tracks cannot be enforced, whereas the linear nature of the 
magnetic anomalies suggests that such continuity exists. In order to understand fully the spatial 
relationships of magnetization, 3D models are necessary. This was an important goal of our 
original proposal that was not attained. The primary responsibility lies with the PI, largely 
because constant supervision was not possible owing to the distance separating the Pi’s 
workplace and the graduate student at the performing institution. Secondary issues were the lack 
of an available student during the first year of this project and the subsequent ability and 
experience of that student in performing geophysical inversions. To reduce commitments to a 
level that could be managed at the Pi’s ability to interact with the performing institution, the 
program there was descoped to a single graduate student. This senior student carried out the 
work described in the attached reprint to this report. During the period of this contract, other 
groups have produced planform-magnetization models {Prurucker et ai, 2000; Arkam-Hamed, 
2002, Cain et ai, 2002) 

5. STRAIN ASSOCIATED WITH MAGNETIC ANOMALIES 

The magnetizatio models can be used as the basis for estimating the amount of extension 
under meridianal stress. We assume that such extension and intrusion occurred while the 
martian geodynamo was active, but that the intrinsic was field was zero at other times. We next 
assume that the maximum observed magnetization intensity corresponds to a model block that is 
entirely composed of intruded, magnetized, igneous rocks. This is an upper limit, as even the 
most magnetized block could still contain unmagnetized material. It is also implicitly assumed 
that the scale of mixing of magnetized and unmagnetized rock is small compared to the block 
width. Under these assumptions, the extension of each block is proportional to the block 
magnetization, and the total extension over the profile can be computed. For the result shown in 
Figure 5, the extension is 1550 km, or a strain of 59% to produce the present-day profile length 
of 4000 km. While this amount of strain is large compared to that expected from vertical 
tectonics alone, it is only moderately larger than that inferred for deep-seated extension of the 
Valles Marineris (10-30%; Anderson and Grimm, 1998) and yet is small compared to the 
essentially unbounded extension that occurs for plate tectonics. The minimum extension 
estimated" from the models above is 830 km, or 26% strain. 




Figure 8. Schematic illustration of assumed relationship between magnetization intensity and crustal strain. 
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6. GEOMETRICAL TESTS 


The paleomagnetic record of Mars can be used in several ways to test the predictions of the 
plate-tectonic and hemispheric-dichotomy hypotheses. The predictive power of the theory of 
plate tectonics lies in the simplicity of its geometric rules for rigid motions on the surface of a 
sphere {McKenzie and Parker, 1967; Morgan, 1968), and this theory can be applied to Mars. All 
plate movements occur along small circles about some apparent rotation pole, and where a plate 
boundary lies in the direction of motion a transform fault results. Therefore the orientation of 
transform faults is the most robust method of determining the geometry of relative rotation 
between two plates. Grimm and Solomon (1989) applied this method to putative fracture zones 
on Venus. If the crustal magnetization reconstructions described above are sufficiently robust 
and strike discontinuities are observed, then we will test to see if they satisfy the geometrical 
requirements of transform faults. 

Ridges and trenches are not geometrically constrained. Where convergent or divergent 
motions" 5 are perpendicular to a plate boundary, however, that boundary will lie along a great 
circle. Regardless of spreading direction, the spreading rate must be proportional to sin(A), 
where A is the angular distance to the rotation pole. This is the second-best method of 
determining plate geometry. The third method, earthquake focal mechanisms, is less reliable and 
is neither available for Mars nor appropriate to paleomotion studies. 



Figure 7 Map of RMS deviations (degrees) for trial southern poles about which magnetic lineations (connected 
diamonds) are tested as small circles. Best fit is 71 S, 8 E (asterisk); proximity to the geographic pole simply 
confirms the sublatitudinal orientation of the lineations. 
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Working directly with the modeled planform of magnetization, we can examine along-strike 
variations in the widths of anomalies which could be a manifestation of a sin(A) spreading-rate 
variation. We can also test whether segments lie along great circles, which would indicate 
normal spreading. Clearly smaller segments can be arbitrarily well fit to a great circle, whereas 
transform faults and/or oblique spreading could render the test invalid for larger segments. As 
the former dominate apparent ridge directions on Earth, arcs that do not fit great circles should 
show transform offsets if plate tectonics is still satisfied (e.g., Southeast-Indian and Pacific- 
Antarctic Ridges on Earth are grossly east-west, but are strongly segmented). 

While there are a number of tests for plate tectonics, the global dichotomy model as presently 
formulated requires only that the crustal magnetizations all follow small circles. For such a test 
(Figure 1), we find an apparent pole at 71 S, 352 W (71 N, 172 W), with an RMS arc deviation 
of 1.2 degrees (Figure 7). The longitude is not particularly meaningful, as it is determined by the 
distribution of the data, but the latitude and goodness-of-fit confirms the visual impression that 
the lineations follow small circles, apparently around some former source beneath the Vastitas 

Borealis. 

The bilateral symmetry of magnetization that was the hallmark of seafloor spreading is not 
immediately obvious in the martian lineations, although the subtle arc shapes to the positive 
anomaly bands 35 S (concave to the south) and 55 S (concave to the north) suggests symmetrical 
spreading about a rotation pole located close by. We have tested the bilateral symmetry of 
magnetization using established statistical techniques ( Grimm and Solomon, 1989). Typical 
correlations -0.65 (Figure 8) do not strongly support the concept of bilateral, plate-like intrusion 
and spreading. Rather, the magnetic anomalies were probably formed by serial, perhaps spatially 
random magmatic intrusions. 



Fieure 8 Bilateral symmetry test for plate tectonics. Each modeled magnetic block, is assumed to be a spreading 
center and the correlation coefficient calculated around it to the maximum possible distance. 
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7. SYNTHESIS: GLOBAL MODELS 

FOR THE ORIGIN OF MAGNETIC LINEATIONS 

We review all alternative hypotheses for the origin of magnetic lineations on Mars and, by 
process of elimination, assess those that remain most likely. 

Plate Tectonics. This is the original hypothesis of Connerney et al. (1999). Arbitrarily large 
(unconstrained) crustal extension and intrusion is allowed. An unmapped location for 
subduction is required to conserve area. There are strong geometrical constraints to motion for 
rigid lithosphere, including the Tit of lineation arcs to great circles and bilateral symmetry of 
lineations. Neither of these constraints is satisfied for the martian magnetization. Finally, a 
near-constant crustal thickness should be produced by decompression melting, but gravity data 
indicates a strong gradient in crustal thickness in the southern highlands ( Zuber et al., 2000). We 
conclude that earth-like plate-tectonics did not form the magnetic lineations. 

Thar sis Tectonics. The load from Tharsis is known to have deformed the crust to great 
distances; could it also have led to the magnitude and direction of the magnetic lineations? Two 
prominent lineations can indeed be fit to great circles with centers near Tharsis, but the anomaly 
pattern as a whole does not follow Tharsis radials. Furthermore, strains from such loading, or 
“vertical tectonics,” are expected to be much smaller, in the range of a few percent; recall that the 
strain associated with the lineations is of order ten percent or more. 

Spin Tectonics. This includes faulting due to the shift of the rotational bulge during spin-axis 
reorientation or polar wander and the relaxation of the rotational bulge due to tidal despinning. 
For the former, we expect N-S oriented extensional structures near the paleopole and only 
moderate strains (several %, e.g., Grimm and Solomon, 1986). For the latter, E-W oriented 
fractures can occur, but they are restricted to poleward of 36° from the equator (. Pechmann and 
Melosh, 1 979). The magnetic lineations of Mars fail the geometrical tests for both of these 

models. 

Northern Crustal Depletion. Subduction or other crustal recycling is one hypothesis for 
formation of the global dichotomy. In this scenario, “back-arc” tensional stresses might be set up 
in the south, which would cause strong gradients in extension, intrusion, magnetization. The 
expected strain magnitude unknown, but it would likely be bigger than vertical or spin tectonics 
but less than plate tectonics: this is consistent with moderately large strains inferred for the 

magnetic lineations. 

Southern Crustal Accretion. This is an alternative hypothesis for formation of the global 
dichotomy in w hich the a megaplume results in generation of a massive crustal plateau. This is 
consistent with the basaltic composition of the southern crust ( Christiansen et al., 1999) and 
indeed is facilitated by a 2.6x gravity-scaled efficiency for decompression melting ( McKenzie , 
1 999). Melt thicknesses of several tens of km result from several tens of km extension (see 
Grimm and Hess, 1997) at hotter Noachian mantle temperatures {Schubert et al., 1992). A 
southern megaplume could generate large, latitudinally varying crustal thickness (80 km at S. 
pole to 35 km at N. plains boundary; Zuber et al., 2000). Radial extension over the uplift causes 


to 


circumferential extension fractures ( Banerdt et al, 1982); subsequent melt intrusion and 
solidification would form the magnetic lineations. This model is most consistent w'ith all 
available data. 

In either the northern-depletion or southern-accretion models, the center of mass (CM) 
becomes displaced from the center of figure (CF), resulting in a spin inbalance. Polar wander 
then brings the CM-CF offset into alignment with the spin axis, i.e., to a N-S configuration. The 
near-perfect N-S alignment of the CM-CF offset has not heretofore been explained. 

Figure 9 (on following page) illustrates the northern-depletion and southern-accretion models 
with subsequent spin-axis realignment. 

8. RELATION OF CRUSTAL MAGNETISM TO SURFACE GEOLOGY 

The spatial distribution of crustal magnetization is statistically correlated to the locations of 
Noachian, southern-highlands valley networks. We hypothesize that at least some of these valley 
networks were carved by hydrothermal discharge over the cooling intrusions represented by the 
magnetization. To test this hypothesis, a numerical model for hydrothermal discharge on Mars 
was formulated and calibrated; a by-product is a better understanding of the physical controls on 
hydrothermal discharge as well as quantitative constraints on crustal permeability imposed by 
geochemistry. We found that the total discharge due to intrusions building that part of the 
southern highlands crust associated with valley networks is comparable to the discharge inferred 
from valley geometry, supporting the hypothesis. Details are given in the attached reprint. 
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[i] Models of hydrothermal groundwater circulation can quantify limits to the role of 
hydrothermal activity in Martian crustal processes. We present here the results of numerical 
simulations of convection in a porous medium due to the presence of a hot intruded magma 
chamber. The parameter space includes magma chamber depth, volume, aspect ratio, and host rock 
permeability and porosity. A primaiy goal of the models is the computation of surface discharge. 
Discharge increases approximately linearly with chamber volume, decreases weakly with depth 
(at low geothermal gradients), and is maximized for equant-shaped chambers. Discharge increases 
linearly with permeability until limited by the energy available from the intrusion. Changes in the 
average porosity are balanced by changes in flow velocity' and therefore have little effect. Water/ 
rock ratios of -0. 1 , obtained by other workers from models based on the mineralogy of the 
Shergotty meteorite, imply minimum permeabilities of 10~ 16 nT during hydrothermal alteration. If 
substantial vapor volumes are required for soil alteration, the permeability must exceed 10 m . 
The principal application of our model is to test the viability of hydrothermal circulation as the 
primary' process responsible for the broad spatial correlation of Martian valley networks with 
magnetic anomalies. For host rock permeabilities as low as 10“ 17 nr and intrusion volumes as low 
as 50 km 3 , the total discharge due to intrusions building that part of the southern highlands crust 
associated with magnetic anomalies spans a comparable range as the inferred discharge from the 
overlying valley networks. INDEX TERMS: 1832 Hydrology: Groundwater transport; 1860 
Hydrologv* Runoff and streamflow; 1545 Geomagnetism and Paleomagnetism: Spatial variations 
(all harmonics and anomalies); 5440 Planetology: Solid Surface Planets: Magnetic fields and 
magnetism; 5114 Physical Properties of Rocks: Permeability' and porosity; KEYWORDS: 
Hydrothermal, groundwater, runoff, crustal magnetism, intrusions 


1. Introduction 

[ 2 ] Hydrothermal circulation is an important part of many 
terrestrial igneous, metamorphic, and sedimentary environments 
and has profound geochemical and biological implications. On the 
Earth it accelerates the cooling of magmatic bodies in systems 
ranging from divergent plate boundaries to individual volcanoes 
and frequently produces discharge in the form of hot springs, 
geysers, and submarine vents. There is also evidence for the past 
existence of hydrothermal systems on Mars (see Farmer [1996] for 
a review). Valley networks are associated with structures that are 
able to force groundwater to the surface [Baker et al , 1992]; these 
include some of the younger volcanoes [Gulick, 1998], ancient 
volcanoes, rifts, and impact craters [Tanaka et al. , 1998]. Many 
more observations have led to groundwater discharge as the 
favored method ot runoff production (see Bciker et al. [1992], 
Baker [2001], and Carr [1996] for reviews) [Malin and Carr , 
1999; Malin and Edgett, 2000]. Hydrothermal circulation may also 
have 5 altered the Martian crust and further produced weathering 
products and soil [Griffith and Shock , 1997; Newsom et al., 1999]. 
The greatest consequence of hydrothermal activity on Mars may be 
its ability to sustain life [Shock, 1996; Farmer, 1996]. 
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[3] We present here numerical models of the thermal convection 
of groundwater in a porous host rock due to the presence of an 
intruded magma chamber. An extensive portion of the available 
parameter space is explored in order to quantify the effects that 
magma chamber volume, depth, aspect ratio, and host rock per- 
meability and porosity have on surface discharge. The results of 
this general approach are then applied to some specific issues. 
First, we use geochemical constraints to bound the permeability of 
the Martian crust. Second, we test the hypothesis that hydrothermal 
circulation can explain the putative correlation observed between 
valley networks and magnetic anomalies [Jakosky and Phillips , 
2001]. We suggest that large amounts of water were circulated 
throughout the southern -high lands crust due to magmatic intrusion 
and that the portion of this water discharged to the surface can 
quantitatively account for the valley networks preserved since the 
end of crustal formation. 

2. Model 

[4] Two-dimensional, axisymmetric representations of hydro- 
thermal circulation in a magma chamber and its host rock were 
modeled with the U.S. Geological Survey (USGS) code HYDRO- 
THERM [Hayba and Ingebritsen, 1994, 1997] and its graphical user 
interface HTpost (P. S. Hsieh, USGS, unpublished material, 2000). 
HYDROTHERM can model temperatures from about 0° to 1200°C 
and pressures from 0.5 to 10,000 bars and keeps track of both liquid 
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and gaseous phases of pure water. It solves mass, momentum, and 
energy balance equations expressed in terms of dependent variables 
pressure and enthalpy. The choice of enthalpy over temperature 
allows the thermodynamic state of the fluid to be specified uniquely 
under two-phase conditions. Viscosity and density for a particular 
temperature and pressure are obtained from a look-up table. 

[ 5 ] The momentum balance equation is Darcy’s law, which for a 
single fluid phase i is 


v, = - ^ (V/> + PjgVz), 
M/ 


( 1 ) 


where y t is the Darcy velocity, k is intrinsic permeability, p, is 
dynamic viscosity, p is pressure, p, is density, g is gravitational 
acceleration, and z is depth. The relative permeability k rc [ quantifies 
the reduction of the flow of phase i due to the presence of the other 
phase. The mass balance (continuity) equation for phase i is 

^ (nS, p,)-i-V • (p,Vj) = 0, (2) 

where n is porosity, t is time, and S, is volumetric saturation (^ wa tcr + 
S s ieam ~ L i.e., the medium is fully saturated). The continuity 
equation for the entire system is the sum of the equations for each 
phase. The energy-balance equation for the entire system is 


( 1 - n) p f .hr + S,p { hi 


fV 




-V-AT w Vr = 0 

( 3 ) 


where h is enthalpy, T is temperature, K m is medium thermal 
conductivity, and subscript r refers to rock matrix properties. 
HYDROTHERM solves the equations by performing Newton- 
Raphson iterations on an equivalent finite difference system (with 
the horizontal dimension expressed in radial coordinates) until mass 
and energy residuals fall below specified maximum values. 

[6] The horizontal extent of the host rock is 20 times that of the 
magma chamber and the vertical extent is 20 km, both of w hich are 
sufficient to accommodate flow from all magma chambers studied. 
The right vertical boundary is not intended to represent the limit of 
a horizontally bounded water source, but to provide enough space 
for a realistic response to a local temperature perturbation. As 
much fluid flows through this boundary as is necessary to balance 
the net flow through the top horizontal boundary. The suitability of 
the chosen horizontal extent was tested by a baseline model 
measuring 40 chamber radii across, which yielded almost the 
identical discharge to the original model. Fluid is allowed to cross 
the upper horizontal and right vertical boundaries, while temper- 
ature and pressure are fixed. Recharge from the surface is typically 
<10% of discharge, indicating that the strength of discharge does 
not depend on infiltration of runoff. Adaptive boundary' conditions 
to prevent infiltration for a dry Mars could therefore be neglected. 
Sufficient recharge or discharge occurs through the right vertical 
boundary to conserve the mass of the system. 

[ 7 ] The initial pressure distribution is hydrostatic, with a surface 
value of 1 bar. A small geothermal gradient (0.5°C/km), applied to 
ensure stable decay of the surface discharge, does not otherwise 
affect the model (the effects of substantial geothermal gradients 
were explored, and the results are described below). The lowest 
temperature HYDROTHERM can handle in this model is 10°C, 
and this was used for the surface boundary condition. Both the 
surface temperature and pressure w'ere adopted for numerical 
convenience, and they neither affect the results nor are intended 
to represent Earth-like climatic conditions. The left vertical and 
lower horizontal boundaries are no-flow', and the temperature and 


pressure are free to vary. Grid spacing is ~ 100 m in both directions 
near the magma chamber and at greater horizontal distances 
increases logarithmically. 

[8] Present Martian conditions may include a permafrost layer 
which must first be melted before discharge is produced. Gulick 
[1998] estimated that the time required to melt a 2 km thick layer 
of permafrost was much less than the lifetime of the hydrothermal 
system, making it unlikely that quantities integrated over the 
lifetime of the system, such as total discharge, should be signifi- 
cantly altered. During the late Noachian, when the hydrothermal 
systems proposed here were active, the permafrost layer was likely 
to be much thinner than 2 km, making melting times even shorter. 
A HYDRO THERM simulation that quantifies the role of ice in our 


models is described below. 

[ 9 ] Our baseline model consists of a 50 km 3 chamber emplaced 
at a depth of 2 km below the surface into host rock of permeability 
10" 16 m 2 . The dimensions of the chamber are expressed in terms of 
its aspect ratio (D/H), defined here as the ratio of diameter to 
height. The baseline chamber has D/H = 2. The volume of the 
chamber is taken after that modeled by Hayba and Ingebritsen 
[1997] and is also of similar size to magma chambers found under 
mid-oceanic ridges [ Burnett et al., 1989]. Note that our baseline 
model is intended only as an example model from which we later 
deviate extensively and does not necessarily represent any sort of 


“ideal” parameter values. 

[ 10 ] The chamber is emplaced instantaneously at a temperature 
of 900°C. Because this approach ignores discharge produced during 
supersolidus cooling and the finite intrusion process, it produces 
relatively conservative results. The magma chamber is initially 
impermeable, but as it cools through a brittle -ductile transition 
(BDT) between 400° and 360°C, it is assumed to fracture and 
become as permeable as the surrounding host rock [Hayba and 
Ingebritsen, 1997]. A semilog form is adopted tor this transition, 
wherein the log of the permeability scales linearly with temperature. 
Note that the “permeability” BDT may differ somewhat from the 
classic “deformational” BDT [e.g., Kohlstedt et a/., 1995]. All 
models use a rock density of 2500 kg m -3 , a thermal conductivity' 
of 2.0 W m~ l K -1 , and a porosity of 1%. Deviations from the 
baseline model include host rock permeabilities of KT 1 and 
10 -15 m 2 , porosity of 25%, magma chamber aspect ratios of 0.2 
and 20, depths of 8.5 and 1 5 km, and volumes of 100 and 2000 km . 

[n] Steam and water fluxes, while fundamental to discharge 
calculations, may also be used to test hypotheses regarding the 
possible geochemical alteration in the system. For each time step, a 
measure of water-to-rock ratio for reactions above a particular 
temperature threshold may be calculated by measuring the total 
mass of water that passes through the region warmer than the 
vuluf* nnrl dividing hv the volume of the region. 


3. Results 

[ 12 ] In the baseline model (magma chamber depth is 2^ km, 
volume is 50 km 3 , DiH = 2, and host rock permeability is 10 m ; 
Figure 1) an initial peak in the surface discharge occurring at only a 
few hundred years following magma emplacement is due to thermal 
pressurization [Delaney, 1982]. Its peak value is not significantly 
greater than discharges that occur later in the model. The extremely 
short-lived nature of this effect, w r hich is less extreme for higher 
host rock permeabilities, contributes negligibly to the total, time- 
integrated mass of water produced by the system (total discharge), 
which we assume to be of primary importance in valley erosion. 
Additionally, thermal pressurization may be weaker in the more 
realistic case of a finite duration intrusion process. 

[!3] The broad peak at 45 kyr is due to thermal convection of 
groundwater. At first, the chamber is supercritical and imperme- 
able, allowing convection in the surrounding host rock only. 
Surface discharge peaks w'hen the magma chamber has cooled 
sufficiently to become permeable and admit adveetion. When only 
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Figure 1. Surficial discharge from the baseline model (magma chamber depth 2 km, volume 50 km , aspect ratio 2, 
host rock permeability 10 _lf m 2 ). The regions separated by solid vertical lines in the enlarged portion of the figure 
denote the different phases present within in the model. The magma chamber becomes permeable over the interval 
0.7-57 kyr. The noise on the curve is due mainly to discretization effects in the presence of steam, which has much 
greater velocity than the liquid phase. 


the liquid phase remains, the exponential decay of discharge 
becomes much smoother and dies out at about 2.8 x 10 years. 

3.1. Permeability 

[h] For simplicity, permeability was homogeneous in all mod- 
els. Gulick [1998] used a value of 10" H nr for hydrothermal 
systems on Mars inferred from young, near-surface Hawaiian 
volcanics. Ingebritsen and Sanford [1998], however, report that 
permeabilities in the east rift zone at Kilauea, while high near the 
surface (1 0' l0 - 1 0 -9 m 2 ), are much lower at depths of just 1 - 2 km 
(I0 -l6 -10“ 15 m 2 ) in rock of the same composition. Manning and 
Ingebritsen [1999] estimate values of between 10“ 17 and 10 4 m 

for the mean continental permeability between 1 and 10 km depth. 
We modeled permeabilities from I0“ 17 to 10 -15 m 2 (Figure 2). The 
upper limit is for computational convenience; we demonstrate 
below that results for higher permeabilities can be extrapolated 
from this range. 

[15] Discharge from the low-permeability (10 m ) model has 
a characteristically large thermal pressurization peak followed by a 
weak main peak. Conduction is the dominant mode of heat transfer 
in this model, and the spatially integrated surface discharge at any 
given time (henceforth called “instantaneous discharge”) is only 
about a tenth that of the baseline model. In the baseline model, 
convection and conduction play comparative roles in heat trans- 
port, while in the high-permeability (10 -15 m 2 ) model, convection 
is dominant. Discharge in all three models lasts between 2500 and 
3000 kyr, but rate of decay of discharge is, in general, proportional 
to permeability, with the 10 17 m 2 model decaying to a tenth of its 
peak value in 200 kyr, the baseline model in 175 kyr, and the 10 
m 2 model in 135 years. 

[16] An important feature of the k - 10 m model is the 
presence of a steam-dominated zone above the magma chamber 
during the first several thousand years. This phenomenon may 
have implications for chemical alteration, as described in the 
discussion below. 


[17] The relationship of discharge to permeability was inves- 
tigated using the high-permeability results of Gulick [1998] as a 
starting point. She modeled a 10 1 m 2 hydrothermal system w r hich 
included the magma chamber implicitly through a heat flux 
boundary condition based on the analytical solution to the con- 
ductive heat flow through the wall of an infinitely long cylinder. 
This implies that no heat was lost through the roof and floor of the 



Figure 2. Surficial discharge from a 2 km deep magma chamber 
of volume 50 km 3 for host rock permeabilities of 10‘ , 10 , 

and 10” 15 m 2 . The total discharge of the k = 10 16 m“ model, in 
kilograms, appears above its curve. The numbers above the other 
curves denote the fraction of this discharge that their correspond- 
ing models produced. 
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Figure 3. Discharge results from HYDROTHERM models of a 
geometry similar to those of Gulick [1998]. An analytic fit to 
these data, with its large permeability asymptote, is shown. The 
discharge obtained by the single permeability modeled by Gulick 
is provided for comparison, as are the total discharges of our 
three full geometry models (from Figure 2). The value obtained 
by Gulick was converted from volume to mass for the purpose of 
this figure, with an assumed water density of 1000 kg m . 


chamber and that no BDT was encountered on cooling. By 
adopting these limitations in HYDROTHERM, the problem is^ 
more efficient numerically and allows us to extend the range of 
assumed permeability (Figure 3) and to compare the results directly 
to those of Gulick [1998]. Our results agreed fairly well, although a 
small difference resulted owing to our more realistic magma 
chamber heat flux and possibly other factors such as numerical 
solution. Overall, we found that the relationship between total 
discharge and permeability could be represented approximately by 

D — a - be~ cK , (4) 


where D and K are the base 10 logarithms of total discharge and 
permeability, respectively, and a, b, and c are positive constants 
(w ith best fit values of 38, 0.089, and 0. 1 2, respectively). As K -> oc, 
D -> a, meaning that as permeability tends to very high values, 
discharge increases asymptotically toward a finite maximum value. 
This is not surprising, since the forces driving the flow are limited 
by the amount of heat contained in the magma chamber. For the 
range of permeabilities covered with our own, more complete 
geometry (K < —15), D is approximately linear in K. Note that the 
total discharges (whose relative magnitudes are given above each 
curve in Figure 2) exhibit this linear behavior. 

[ig] For our range of permeabilities, total discharge is approxi- 
mately the same for both geometries. At higher magma chamber 
aspect ratios, however, it is no longer reasonable to assume that 
heat is lost through the chamber walls alone, and the total 
discharges for the two geometries are expected to diverge. Instan- 
taneous discharge is not the same for both geometries, even with 
DiH - 2. Our geometry produces a greater maximum than the 
simpler geometry (~ 1 .75 times as high iork= 1 0“ 16 m 2 ) but drops 
more rapidly thereafter. 

[ 19 ] Magma chamber cooling times, defined here as the average 
time taken for chamber nodes to cool below a specified threshold 
temperature, are of interest in these models. For all three perme- 
abilities the cooling time for a threshold of 450°C (i.e., half the 
emplacement temperature) is ~30 kyT. While an increase in host 
rock permeability increases the velocity of the flow alongside the 


still impermeable magma chamber, it does not significantly 
enhance cooling [see Norton and Knight , 1977]. Once the chamber 
becomes permeable, however, cooling progresses as different rates 
in each model. The time taken for the chamber to cool to 250°C 
was 67, 61, and 40 kyrs for permeabilities of 10” 1? , 10 16 , and 
10” 15 m 2 , respectively. 


3.2. Volume 

[20] Head and Wilson [1994] suggest that Martian magma 
reservoir volumes could be as great as 2000 km 3 and that chamber 
depths are most likely to range from 8 to 12 km. We ran models of 
DiH =2 magma chambers with volumes of 100 and 2000 km 3 at a 
depth of 8.5 km and with host rock permeability of 10” rrr 
(Figure 4). The 100 km 3 chamber produced ~2.5 times as much 
discharge as the 50 km 3 chamber, while the jump from 100 to 2000 
km 3 resulted in a factor of 73 increase. This approximately linear 
relationship is reflected in the magnitude of the instantaneous 
discharge, whose peak value has proportional increases. This, 
coupled with an increase in cooling time with volume (and there- 
fore an increase in the discharge lifetime), explains the observed 
total discharge increase. 

[ 21 ] Although not obvious on the semilog plot of Figure 4, the 
discharge produced by the 2000 km 3 chamber drops off at ■~800 
kyr, which is in agreement with values obtained by Cathles et al. 
[1997] for a 2500 km 3 chamber with host rock permeabilities 
ranging from 4 x 10” 17 to 10” 16 m 2 . 


3.3. Depth 

[ 22 ] Models with 50 km 3 magma chambers at depths of 8.5 and 
15 km and with host rock permeability of 10” m were run 
(Figure 5). Total discharge decreases by a factor of ~ 1.5 with each 
6.5 km increase in depth. Similar depth dependence exists for the 
same three chamber depths at host rock permeabilities of 1 0 and 
!0~ ]5 m 2 . A summary of the total discharge of all nine depth and 
permeability combinations is depicted in Figure 6. The very small 
geothermal gradient in these models makes magma chamber cool- 
ing times only marginally sensitive to changes in depth. All three 
chambers cool to half their emplacement temperature (450°C) in 
~30 kyr. The 2 km deep chamber cools to 250°C in 61 kyr, while 
1 c tnirpc 71 1 a/r Flow from deener chambers must 



Figure 4. Surficial discharge from magma chambers of volume 
50, 100, and 2000 km 3 , all at a depth of 8.5 km, an aspect ratio of 
2, and with host rock permeability of 10” 16 m 2 . The total discharge 
of the 50 km 3 model appears above its curve. The numbers above 
the other curves denote the fraction of this discharge that their 
corresponding models produced. 
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Figure 5. Surficial discharge from 50 km 3 magma chambers at 
depths of 2, 8.5, and 15 km, all with aspect ratio 2, and with host 
rock permeability ot 10 1<: m~. The total discharge of the 2 km 
deep chamber appears above its curve. The numbers above the 
other curves denote the fraction of this discharge that their 
corresponding models produced. 


Figure 7. Surficial discharge from 50 km 3 magma chambers with 
DIM = 0.2, 2, and 20, all at a depth of 2 km, and with host rock 
permeability of 10 16 m 2 . The total discharge at DIH = 2 appears 
above its curve. The numbers above the other curves denote the 
fraction of this discharge that their corresponding models 
produced. 


travel farther in order to expel heat from the system and is therefore 
less efficient at cooling the magma chamber. 

3.4. Aspect Ratio 

[23] Magma chambers with DIH = 0.2, 2, and 20 were modeled 
(Figure 7). These correspond to chamber radii of 1.2, 2.5, and 5.4 
km and chamber heights of 12, 2.5, and 0.54 km. Chamber depth 
in all cases is 2 km, and chamber volume is 50 km . The horizontal 
extent was fixed at 50 km for all three models, rather than scaled 
with the chamber radius, so that the influence of the right vertical 
boundary was the same in each model. The total discharge of these 
models is controlled by the cooling time of the magma chamber, 
and therefore its surface area/volume ratio AIV (1.89, 1.59, and 
4,06 km -1 for the three aspect ratios, respectively). There is an 
approximately linear inverse relationship between A V and total 
discharge; that is, high AIV implies low total discharge and vice 
versa. The effects of aspect ratio are observed more directly when a 
substantial geothermal gradient is applied (see section 3.5). 



Figure 6. Total discharge from 50 km 3 magma chambers of 
aspect ratio 2 at depths of 2, 8.5, and^!5 km and host rock 
permeabilities of 10 17 , 10 1C , and 10 m . 


[24] In all three models the flow pattern consists of a single, 
clockwise rotating convection cell alongside the chamber. For 
DIH = 20 this pattern does not give way to a series of cells above 
the chamber roof, as one might expect for flow between two 
infinite horizontal surfaces at different temperatures. 

3.5. Geothermal Gradient 

[25] Models with magma chamber aspect ratios of 0.2, 2, and 20 
were run with an initial host rock geothermal gradient of 20°C/km, 
perhaps representative of early Mars [ Schubert et ai, 1992]. Die 
presence of a geothermal gradient significantly affects the hydro- 
thermal discharge. The flat, sill-like chamber (. DIH - 20) produces 
the greatest peak discharge because of its large horizontal expo- 
sure. The discharge dissipates more rapidly, however, because the 
chamber, having the greatest AIV and being oriented perpendicular 
to the main flow direction, cools rapidly. Conversely, the tall, pipe- 
like chamber (DIH = 0.2), being immersed in warmer temperatures 
at depth, cools more gradually. Its geometry produces the largest 
convection cell and offers minimum obstruction to flow, resulting 
in the greatest total discharge despite its low peak value. Overall, 
tenfold variations in aspect ratio produce changes in discharge of 
less than a factor of 3. 

[26] The effect of geothermal gradient is also observed in cool- 
ing times. The times taken for magma chambers to cool to half their 
emplacement temperature are, in order of increasing aspect ratio, 
25, 45, and 4.0 kyr, respectively. In the absence of a significant 
geothermal gradient, the same chambers cool in 16, 28, and 3.5 kyr, 
respectively. An 8.5 km deep chamber in a geothermal gradient of 
20° C/km cools to half its initial temperature in 52 kyr (as opposed to 
30 kyr for the model with negligible geotherm). It also produces ~7 
times as much discharge (with a peak 3.5 times higher) as the 
identical model with a negligibly small geothermal gradient. ^ 

[27] It should be noted that for permeabilities >10 1 m , the 
Rayleigh number of a plane porous medium [e.g., Turcotte and 
Schubert , 1982] at 20°C/km indicates that weak convection may 
occur in the absence of a magma chamber. This phenomenon is 
observed in our high geothermal gradient models when no magma 
chamber is emplaced. Free convection in the terrestrial crust has 
been invoked by Raffensperger and Garven [1995] to explain the 
location of uranium ore deposits in sedimentary basins in Canada 
and Australia. Travis el al. [2001] showed that free convection may 
be capable of melting significant volumes of subsurface ice in the 
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Figure 8. Baseline model with an initially 1 km thick subsurface ice layer shown at (left) 15 kyr and (right) 59 kyr. 
Both horizontal and vertical dimensions in each image are approximately 5 km. Water and steam velocities are shown 
by black and gray streamers, respectively. Only every other row and column are shown. Shaded regions (apart from 
the ice) represent different fluid phases as described in the key. At 59 kyr, only a small remnant of steam remains, and 
the magma chamber has become as permeable as the host rock. 


Martian crust. However, all of these models, including our own, 
have homogeneous permeability in the converting zone: realistic 
vertical heterogeneity in planetary crusts will both inhibit the 
development of large-scale crustal convection and decrease the 
efficiency of heat transfer. Although background geothermal gra- 
dients may affect cooling times of intrusions, we view their 
contribution to hydrothermal circulation as doubtful. 

3.6. Ice 

[28] Gulick [1998] suggested that a 2 km thick subsurface 
permafrost layer above a 50 km 3 magma chamber would melt in 
a few tens of thousands of years. This is of the same order as the 
time taken for discharge to peak in our own models, indicating that 
ice could significantly reduce total discharge. We thus ran the 
baseline model with a 1 km thick layer of subsurface ice (Figure 8), 
modeled after Bonacina et ai [1973], who approximated the 
melting process as a cooling period over a small finite temperature 
interval AT During melting the material is assigned an augmented 
specific heat given by 

Ci = Ci+ Sr’ 

w r hcre subscripts S and L refer to the solid and liquid phases, 
respectively, and L is the latent heat of fusion. We used c L = 1000 J 
kg -1 K' 1 , 1 = 3.34 x 10 5 J kg' 1 , and r=5°C, giving c L = 6.7% x 
10 4 J kg” 1 K" 1 . 

[29] The hydrothermal system took 52 kyr to melt a hole in the 
ice. The strong upward convection associated with increasing 
magma chamber permeability did not noticeably increasing melt- 
ing rate. The total discharge produced by the model was 2.95 x 
10 u kg, about a quarter that of the baseline model. This relatively 
severe reduction (which is expected to be worse for deeper 
chambers and lower host rock permeabilities) places an upper 
bound on permafrost thickness during valley formation on Mars. 
For valleys to form through sapping processes alone, the perma- 
frost must be melted through, or aquifers carrying groundwater 
beneath the permafrost must intersect the surface. In either case, 
the permafrost can be no thicker than a few hundred meters. 


3.7. Water Table 

[ 30 ] Surface discharge will occur only if convection is strong 
enough to elevate groundwater from the initial water table depth to 
the surface. We compared with hydrostatic conditions the vertically 
integrated product of density, gravitational acceleration, and depth 
in order to estimate the height a water table may attain through 
thermal-convective expansion. Integrations were performed over 
all vertical columns of finite difference blocks at all times and for a 
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Table 1. Water/Rock Ratios by Mass Calculated for 8.5 km Deep Chambers 


of the Indicated Volumes and Host Rock Permeabilities k 


Threshold 

/fc= 10“ 15 m 2 

k = 10" 16 ra 2 

k - 10~ 16 m 2 

Temperature, °C 

Volume = 50 km 3 

Volume = 50 km 3 

Volume = 2000 km 3 

150 

0.12 

0.025 

0.067 

200 

0.18 

0.037 

0.10 

250 

0.28 

0.054 

0.15 


range of water table depths. The discharge was calculated from 
those parts of the water table that intersected the surface and was 
compared to zero depth water table models. The results (Figure 9) 
show that models with a water table 50 m deep produce only 10% 
of their zero depth water table equivalents, while models with a 
1 00 m deep table produce only 2%. 

4. Discussion 

4.1. Implications for Hydrothermal Alteration 

[31] Water/rock ratio ( W/R ) can be an important influence on 
the mineralogy of alteration products in hydrothermal systems. 
Martian hydrothermal alteration is thought to occur at low W/R 
(<10 by mass) [ Griffith and Shock , 1997; Newsom et ai, 1999], in 
which case its effect is small, and the initial mineral composition is 
the primary influence on the alteration assemblage [Griffith and 
Shock , 1997]. 

[ 32 ] HYDROTHERM mass flux results allow the W!R of a 
hydrothermal model to be calculated. At each time step, regions of 
the model above a specific threshold temperature are identified, 
and the total flux passing through them is calculated. On the basis 
of these data, a W/R (by mass) is calculated for each finite 
difference block in the model that is at some time above the 
threshold temperature. The average W/R values for the three 
reaction temperatures modeled by Griffith and Shock [1997], and 
for various model dimensions, are shown in Table t. 

[ 33 ] The large velocities in the k — 10 m model give the 
highest W/R for all threshold temperatures, ^5 times that of the k - 
10“ 16 m 2 model with a chamber of the same volume. Negligible 
geothermal gradient and low host rock permeability in our models 
are both factors that contribute to small W/R . Typical Martian 


geotherms (>50°C/km) will lead to larger W/R, so these estimates 
can be viewed as lower bounds. The presence of a significant 
geotherm may be expected to increase discharge by a factor of 
about 3 or 4 (as observed in the model with an 8.5 km deep 
chamber with geothermal gradient), and since W/R is directly 
proportional to the discharge flowing through the alteration area, 
it may experience a similar increase. This would still leave most 
values listed in Table 1 below unity'. 

[ 34 ] A primary goal of the work by Griffith and Shock [1997] 
was to calculate the amount of water bound to the rock during 
alteration. In a model based on the Shergotty SNC meteorite, ~8 A o 
of the final mineral assemblage (by weight) was water, implying 
that water/rock ratios >0.08 would be necessary to sustain hydro- 
thermal circulation. Assuming the Shergotty composition is suffi- 
cientlv generic to be applied to our own models, a ff/ R limit of 
0.08 implies a minimum permeability of 10 1 nT for a 150 C 
reaction in host rock surrounding a 2000 km 3 chamber (Table 1). 

[ 35 ] Newsom et ai [1999] suggest that the relative abundances 
of mobile elements such as sulfur, chlorine, sodium, and potassium 
may be explained by the presence of a mixture of neutral -chloride 
and acid-sulfate fluids during soil formation. Production of the 
latter fluid requires vapor transport [Rye et ai, 1992], and Inge- 
britsen and Sorey [1988] discuss situations in which vapor-domi- 
nated zones may occur. Their models require combinations of low- 
permeability barriers and in some instances topographic gradients 
to sustain vapor-dominated zones, which develop in the shallow 
subsurface only. These specialized structures have not been 
included in our generalized models; nonetheless, our highest 
permeability model (1CT 15 m 2 ) with the shallowest chamber (2 
km) does produce a short-lived (few thousand years) two-phase 
zone between magma chamber and surface. Steam develops here 
because of a combination of low pressures (which drop to a 



0 150 300 450 600 nT 


Figure 10. Overlay of the absolute values ofthe 200 km vertical magnetic anomalies [Purucker etal 2000] and the 
valley networks [Kieffer, 1981] (cylindrical projection). Valleys were not mapped for latitudes below- 60 S. e 
anomalies range from 0 to 670 nT. 
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Table 2. Statistical Data for the Correlation Between Valley Networks and Magnetic Anomalies 


Magnetic 






Threshold, nT 

Po 

P 

(p - PqVp 

d m = 2° 

d„, - 3 

] 

0.896 

0.962 

0.069 

0.0785 

0.182 

2 

0.827 

0.928 

0.109 

0.0204 

0.102 

5 

0.669 

0.827 

0.191 

0.00228 

0.0480 

10 

0.469 

0.658 

0.287 

0.00223 

0.0345 

20 

0.295 

0.465 

0.366 

0.00459 

0.0590 

50 

0.126 

0.201 

0.373 

0.0686 

0.303 

100 

0.0521 

0.0818 

0.363 

0.260 

0.450 


0 


d„ = 4° 

0.494 

0.210 

0.172 

0.160 

0.166 

0.330 

0.575 


d m = 6° 

0.464 

0.654 

0.268 

0.434 

0.343 

0.612 

1 


minimum of about two thirds hydrostatic pressure) and high 
temperatures. At lower permeabilities, fluid rising from the cham- 
ber does not transport enough heat upward to elevate temperatures 
to the required levels for steam production. 

[ 36 ] If substantial quantities of vapor are required for soil 
alteration, then the abundance of steam in our high-permeability 
model may point to 10“ 15 m 2 as a lower bound to permeability, a 
somewhat tighter constraint than that imposed by WIR ratio alone. 
Alternatively, if volcanic aerosols [Newsom et al. , 1999] or basalt 
palagonitization [McSween and Keil , 2000] produce sufficient soil 
alteration, then no such limit is necessary. 

4,2. Relationship of Hydrothermal Circulation 
to Valley Networks and Magnetic Anomalies 

[ 37 ] The apparent spatial correlation between valley networks 
and magnetic anomalies suggests a link between processes involv- 
ing the acquisition of thermoremnance in a cooling intrusion and 
processes involving the production of surface water for the creation 
of fluvial channels, i.e., hydrothennal circulation. An overlay of 
maps (Figure 10) showing the valley networks [Kieffer, 1981] and 
the vertical component of the magnetic anomalies as measured by 
the Mars Global Surveyor magnetometer [Acuna et al. , 1999] 
inferred at a constant altitude of ~200 km [Purucker et al., 
2000] visually suggests that the two are correlated [ Jakosky and 
Phillips, 2001]. 

[ 38 ] A statistical analysis of the relative positions of valley 
networks and magnetic anomalies quantifies the correlation. A 
binomial test was performed, with success defined as the occur- 
rence of a valley network and a magnetic anomaly (w ith vertical 
component above a specified threshold) in the same data bin. We 
consider only the southern highlands in these calculations. The 
valley network map [Kieffer, 1981] contains the coordinates of the 
0.25 by 0.25 degree bins that contain valleys; the number and 
length of the valleys are not quantified. The magnetic anomaly map 
[Purucker et al., 2000] is binned at 1 by l degree and was 
regridded at 0.25 degree to match the valley grid. The number of 
degrees of freedom must, however, be adjusted to match the true 
spatial resolution, which is limited by the magnetic data. The total 
number of 0.25 degree bins is therefore multiplied by the quantity 



where d m is the length scale in degrees for magnetic resolution. We 
determined that the variogram for the 200 km magnetic field 
reaches half of its sill (asymptotic) value in 200 km and 90% of the 
sill value in 400 km, in agreement with the rule of thumb for 
potential fields that the spatial resolution is approximately equal to 
measurement altitude. Therefore d m ~ 200 km (-3 degrees) may be 
most appropriate, but we consider a range from 2 to 6 degrees. 

[39] Relevant probabilities may be obtained using the total 
number N of 0.25 by 0.25 degree bins in the region of interest, 
the number n of these bins containing valleys, the number m of 
bins with magnetic anomalies above a specified threshold, and the 


number C of correlations. The observed correlation probability is 
then p = Cln, while the expected correlation probability for a valley 
placed randomly in the region of interest is po = m/N. The 
probability q of obtaining C or more correlations if the n valleys 
are placed randomly in the region of interest can be obtained using 
the cumulative binomial distribution with zC, zn, and pQ as the 
minimum number of successes, the number of independent trials, 
and the probability of a single success, respectively. A low value of 
q implies that the relative distribution of valleys and magnetic 
anomalies on Mars is not what one would expect from chance. 

[ 40 ] The results for a range of magnetic anomaly thresholds 
(Table 2) indicate that minimum chance probabilities occur near a 
10 nT threshold for the vertical magnetic field at 200 km altitude. 
The area containing magnetic anomalies tails off sharply at higher 
thresholds, allowing a greater probability of chance correlation. At 
lower thresholds the entire map is considered anomalous, and so 
again it is easier to produce randomly the observed correlation. 
Over length scales at which the magnetic anomalies can be 
described as strongly coherent {d m < 4 degrees), the probability 
of a chance correlation is relatively small (q < 0.16). 

[41 ] A genetic correlation requires that the valley networks and 
magnetic anomalies are the same age. A significant majority of 
valleys are Noachian (70-92% [Scott and Dohm , 1992; Carr , 
1996]); many of those that are younger are not in the southern 
highlands and so are excluded from our survey. The magnetic 
anomalies, because of their inferred deep-crustal origin (see prop- 
osition 1 below) without surface manifestation, must also be very 
old. However, even the oldest valley networks individually pre- 
serve only some part of the Noachian that was not subsequently 
locally erased, whereas the magnetic anomalies probably reflect a 
greater span of crustal history. In other words, some valley net- 
works that were associated with magnetic anomalies may have 
been resurfaced, whereas other valley networks may have formed 
subsequent to emplacement of the magnetic anomalies by hydro- 
thennal or other processes. The normalized excess of correlated 
valleys and magnetic anomalies {p - po)/po (Table 2) may be taken 
as representative of the fraction of valleys for which a genetic 
correlation may be inferred, between about one quarter and one 
third. 

[ 42 ] The role of hydrothermal circulation in the relationship 
between valley networks and magnetic anomalies may now be 
described in terms of a central hypothesis composed of two main 
propositions, defined and discussed in the following sections. 

4.2.1. Proposition 1. [ 43 ] The first proposition is that the 
magnetic anomalies formed as intruded crust and that the 
acquisition of thermal remnant magnetization (TRM) occurred at 
relatively great depth. The strong observed magnetizations in the 
Martian cmstof20-40 Am 2 [Connerneyetal, 1999; Grimm, 2000] 
imply magnetization depths of up to a few tens of kilometers. The 
presence of magnetic anomalies in Arabia Terra, which has been 
strongly resurfaced [ McGill , 2000; Hynek and Phillips , 2001], also 
points to a deep origin. On Earth, Layer 3 gabbros of the oceanic 
crust are magnetized to the extent that they contribute between 25 
and 75% of the observed marine anomalies [Pariso and Johnson, 
1993], while the underlying mantle is unmagnetized [Wasilewski et 
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al . , 1979], Therefore deep crustal magnetization of Mars is 
reasonable. 

[ 44 ] The mineral composition of the magnetized material pro- 
ducing the magnetic anomalies is likely to contain magnetite or 
hematite as the primary magnetic carrier. Although magnetite is 
generally favored, there is much support for hematite [e.g., Con- 
nerney el al., 1999]. Kletetschka et al. [2000] show that for an 
applied magnetic field of 0. 1 mT (about twice the strength of the 
Earth’s present geomagnetic field), multidomain hematite reaches 
maximum TRM saturation, whereas magnetite reaches only a few 
percent thereof. 

4 . 2 . 2 . Proposition 2 . [ 45 ] The second proposition is that 
hydrothermal discharge attending crustal formation processes in 
the southern hemisphere of Mars was sufficient to provide the water 
necessary to carve the valley networks. We assume that where 
instantaneous discharges predicted by our models are too small to 
do significant geomorphic work, topographical variations and near- 
surface heterogeneity' in the host rock permeability (especially in the 
horizontal dimensions) are sufficient to concentrate discharge to the 
required levels. Hydrothermal systems on Earth exhibit such 
discretized concentration of outflow, as evidenced by the presence 
of geysers and springs rather than diffuse outflow everywhere 
above the magmatic intrusion. The Martian valley networks are 
characterized by low drainage densities, implying again that crustal 
heterogeneities may localize discharge. 

[ 46 ] Testing this proposition requires knowledge of the amount 
of water necessary to erode the valley networks, the total amount of 
water available to hydrothermal systems, and the actual hydro- 
thermal discharge produced. An estimate of the required water 
volume for valley erosion may be made from values of areal 
coverage, drainage density, valley cross section, and sediment-to- 
water ratio. Using a map of drainage densities [ Carr and Chuang, 
1997], we estimate the total area covered by valley networks to be 
about 1.4 x 10 7 km 2 . Since we estimated earlier that only one 
quarter of the valley networks may presence direct interaction with 
hydrothermal systems, we use a reduced effective area of 3.6 x 10 
km 2 . Carr and Chuang [1997] calculated a globally averaged 
drainage density of 0.0032 km' l . The product of effective area and 
drainage density, multiplied by typical valley width and height (5 
km and 150 m, respectively), results in 8.6 x 10 3 km 3 of removed 
material. Sediment-to-wuter ratios may range from 1:4 to 1:1000 
[ Gulick , 1998, and references therein], implying that the volume of 
water required to erode the valleys was between 3.5 x 10 4 and 8.6 
x 10 6 km 3 , equivalent to a global water layer between 0.2 and 60 
m deep. These values are well below the estimate of hundreds of 
meters for the global crustal inventory, indicating that discharge 
from the valleys had a negligible impact on the global water 
budget. They further imply that discharged water need not have 
been recharged to the crust. 

[ 47 ] The total hydrothermal discharge produced is computed as 
follows. First, we assume that each magma chamber contributing 
to crustal fonnation intrudes into steady ambient temperature and 
pressure conditions; this is most likely if intrusions that formed the 
southern highlands moved between different loci rather than 
spreading from a single location [Grimm, 2000]. In this way, the 
discharge contribution from a single intrusion is just that of its 
equivalent HYDROTHERM model. We assume further that the 
crust covering the area occupied by valley networks (as calculated 
above) is eventually built up to a depth of 20 km by magma 
chambers packed side by side and one on top of the other. We sum 
the total discharges from each intrusion to calculate the total mass 
of surface water produced. We consider only those intrusions that 
contribute to magnetic anomalies, discarding other intrusive 
events. Determining the relative contribution of individual intru- 
sions is not possible, but the probability p 0 may be used as an 
indicator of the appropriate fraction to be considered. Its value for a 
magnetic threshold of 10 nT, i.e., 0.469 (Table 2), is applied to all 
of our results. 
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Figure 11 . Summary of global discharge production (by mass, in 
kg) for various model parameters. The two dashed lines indicate 
the bounds on the mass of water required to erode the valley 
networks. Represented are discharge values for different magma 
chamber aspect ratios, for host rock permeabilities, for magma 
chamber volumes, and for the baseline model with different 
assumed water table depths. Note that the baseline model {DIH ~ 2, 
k = 10' 16 m 2 , chamber volume = 50 km 3 , water table depth = 0 m) 
appears in each vertical group of bars. 

[48] A summary of discharge production over the entire region 
of interest is shown for different models in Figure 11. The two 
thick dashed lines denote the bounds on total required discharge. 
Represented are discharge values for different magma chamber 
aspect ratios (all with a volume of 50 km 3 , host rock permeability 
of 10“ 16 m 2 ), host rock permeabilities (all chambers have DtH- 2 
and a volume of 50 km 3 ), magma chamber volumes (all chambers 
have D/H = 2; host rock permeability is 10“ 16 m 2 ), and values for 
the baseline model with different assumed water table depths. 

[ 49 ] All models meet or exceed the discharge production 
requirements, although it should be noted that if the k — 10 
m 2 model is assumed to have a water table of 100 m, its production 
will fall below the minimum required value. Other factors such as 
evaporation may further reduce the discharge available to carve 
valleys. The present evaporation mass flux is likely in the region of 
4 x 10“ 8 kg m -2 s' 1 [Ingersoll, 1974], which would reduce the 
effective discharge by as much as an order of magnitude or more. 
Ice may also reduce discharge (section 3.6). A 1 km thick layer of 
subsurface ice in our baseline model causes a 75% drop in total 
discharge. In models with other parameter values (e.g., greater 
magma chamber depth and smaller host rock permeability), ice 
may inhibit discharge more severely, if not completely. 

5. Conclusions 

[ 50 ] In numerical models of Martian hydrothermal systems we 
explored the confrol on surface discharge of magma chamber 
depth, volume, aspect ratio, and host rock permeability and 
porosity. Discharge has an approximately linear relationship to 
magma chamber volume and host rock permeability (in the range 
1 0“ 17 — 1 0 — 1 5 m 2 ). The influences of depth and aspect ratio are 
weaker, and that of porosity is negligible. 

[ 51 ] Some geochemical aspects of Martian hydrothermal sys- 
tems were considered by calculating water/rock ratios in our 
numerical models at various reaction temperatures. Ratios^tend to 
be low' but sufficiently large at mean permeabilities >10' 1 m for 
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groundwater flow to be sustained, consistent with the expected 
storage of water in alteration assemblages. The presence of a short- 
lived vapor-dominated zone in our model with high host rock 
permeability (10“ 15 m 2 ) and shallow chamber depth (2 km) 
suggests that hydrothermal alteration processes may be responsible 
for the observed relative abundances of certain salts in the Martian 
soil, although other forms of alteration should not be excluded. 

[52] Crustal formation processes which formed the magnetic 
anomalies observed on Mars today may have been attended by 
hydrothermal circulation that also provided surface water for valley 
network erosion. This idea is in agreement with the observed 
spatial correlation between magnetic anomalies and valleys and 
was tested further within the framework of a central hypothesis 
made up of two propositions. The first is that the magnetic 
anomalies formed as intruded crust and that the acquisition of 
thermoremnance occurred at relatively great depth. The second is 
that hydrothermal discharge attending global cmstal formation 
processes is sufficient to provide the water necessary to carve the 
planet’s valley networks. 

[53] We tested this hypothesis using the numerical models 
described above, assuming that the crust was formed by the 
heterogeneously spaced and timed intrusion of multiple magma 
chambers, each of which produced the discharge predicted by its 
individual numerical model. We determined that many model 
configurations in the explored portion of the parameter space were 
capable of producing sufficient water to erode those valley net- 
works statisticallv related to hydrothermal circulation. In particular, 
modest cmstal permeabilities of 1(T 16 - 10“ 15 m 2 can produce the 
discharge required to carve valleys and satisfy geochemical con- 
straints, even in the presence of mitigating factors such as evap- 
oration and finite water table depth. 
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